/*	Interactive Gravitational Simulator
	Project: IGS Core
	File: Graviton.h
	Author: Mike Bantegui <mbante2@gmail.com>, Hofstra University 
	Copyright (C) 2012 - 2013 Mike Bantegui

	This file is part of the IGS software package for simulating
	the N-Body problem in real-time. This was developed as part
	of a departmental honors thesis project at Hofstra University
	in the Spring 2012 semester.

	IGS is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation; either version 2 of the License, or
	(at your option) any later version.

	IGS is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/

#ifndef IGS_GRAVITON_H
#define IGS_GRAVITON_H

#include "WorldPoint.h"

namespace IGS
{

struct Body;
struct Node;
struct ViewNode;

// Represents anything which emits or receives the gravitational force
struct Graviton : public WorldPoint
{
	// Mass of the Graviton
	double Mass;
	// Potential evaluated at the Graviton's position
	double Potential;
	// The radius of the Graviton
	double Radius;
	// The acceleration and higher order derivatives of the Graviton
	Derivatives Derivative;

	Graviton();
	Graviton(const double mass, const Vector &position, const Vector &velocity);
	virtual ~Graviton() = 0;

	virtual ViewNode *AsReadOnly() const = 0;

	virtual void CheckTimescale(const double timescaleSquared) = 0;

	virtual void Pairwise(Body &other, const double softening) = 0;
	virtual void Direct(Body &other, const double softening) const = 0;

	// Multipole Expansion
	virtual double MultipolePotential(const Vector &dr, const double r1, const double r2) const = 0;
	virtual void AddMultipole(Node &branch) const = 0;
};

template <class T1, class T2>
void Pairwise(T1 &body1, T2 &body2, const double softening, const Vector &dr, double r2)
{
#ifdef PlummerSoftening
	if (softening > 0)
		r2 += softening;

	const double r1 = sqrt(r2);
#else
	double r1 = sqrt(r2);

	if (softening > 0)
	{
		r1 += softening;
		r2 = r1 * r1;
	}
#endif

	// Potential
	const double p1 = body2.Mass / r1;
	const double p2 = body1.Mass / r1;
	body1.Potential -= p1;
	body2.Potential -= p2;

	Derivatives d1, d2;

	// Acceleration
	const double da1 = +p1 / r2;
	const double da2 = -p2 / r2;
	d1.Acceleration = da1 * dr;
	d2.Acceleration = da2 * dr;

	const Vector dv = body2.Velocity - body1.Velocity;
#if DerivativeCount > 1
	// Jerk
	const double alpha = Dot(dr, dv) / r2;
	d1.Jerk = da1 * dv - (3 * alpha) * d1.Acceleration;
	d2.Jerk = da2 * dv - (3 * alpha) * d2.Acceleration;
#endif

	const double v2 = Dot(dv, dv);
#if DerivativeCount > 2
	// Snap
	const Vector da = body2.Acceleration - body1.Acceleration;
	const double beta = (v2 + Dot(dr, da)) / r2 + alpha * alpha;
	d1.Snap = da1 * da - (3 * beta) * d1.Acceleration - (6 * alpha) * d1.Jerk;
	d2.Snap = da2 * da - (3 * beta) * d2.Acceleration - (6 * alpha) * d2.Jerk;
#endif

	// Accumulate derivatives
	body1.Derivative += d1;
	body2.Derivative += d2;

	const double ctsq = r2 / std::max(std::max(p1 + p2, v2), 1E-4);
	body1.CheckTimescale(ctsq);
	body2.CheckTimescale(ctsq);
}

template <class T1, class T2>
void Direct(T1 &body1, const T2 &body2, const double softening, const Vector &dr, double r2)
{
#ifdef PlummerSoftening
	if (softening > 0)
		r2 += softening;

	const double r1 = sqrt(r2);
#else
	double r1 = sqrt(r2);

	if (softening > 0)
	{
		r1 += softening;
		r2 = r1 * r1;
	}
#endif

	// Potential
	const double p1 = body2.Mass / r1;
	const double p2 = body1.Mass / r1;
	body1.Potential += -p1 + body2.MultipolePotential(dr, r1, r2);

	Derivatives d1;

	// Acceleration
	const double da1 = +p1 / r2;
	d1.Acceleration = da1 * dr;

	const Vector dv = body2.Velocity - body1.Velocity;
#if DerivativeCount > 1
	// Jerk
	const double alpha = Dot(dr, dv) / r2;
	d1.Jerk = da1 * dv - (3 * alpha) * d1.Acceleration;
#endif

	const double v2 = Dot(dv, dv);
#if DerivativeCount > 2
	// Snap
	const Vector da = body2.Acceleration - body1.Acceleration;
	const double beta = (v2 + Dot(dr, da)) / r2 + alpha * alpha;
	d1.Snap = da1 * da - (3 * beta) * d1.Acceleration - (6 * alpha) * d1.Jerk;
#endif

	// Accumulate derivatives
	body1.Derivative += d1;

	// We use square velocity scale that is the largest of:
	// 1) The square of relative velocity between each body
	// 2) The sum of the magnitude of the acceleration of each body
	// 3) A constant factor of 1E-4
	
	// Note: It may also be interesting to observe the dynamics if we consider
	// a time scale factor of |dr|^2 / Dot(v_i, v_j). Bodies with orthogonal
	// velocity vectors will have larger time scale parameters than bodies
	// which have parallel velocity vectors

	const double ctsq = r2 / std::max(std::max(p1 + p2, v2), 1E-4);
	body1.CheckTimescale(ctsq);
}

}

#endif